This document shows how to calculate candidate Pillars and their constituent Indicators at Census Block Group, Tract, and Service Area/ Municipality level. Here, we load the data we pre-processed in the previous exercise
b <- read_sf("../data/boundary.geojson")
bg.l <- read_sf("../data/blockgroups_long.geojson")
bg.w <- read_sf("../data/blockgroups_wide.geojson")
tr.l <- read_sf("../data/tracts_long.geojson")
tr.w <- read_sf("../data/tracts_wide.geojson")
pl.w <- read_sf("../data/place_wide.geojson")
pl.l <- read_sf("../data/place_long.geojson")
lead <- read_sf("../data/parcels_lead_service_lines.geojson")
lead.c <- sf::st_centroid(lead)
calculate_volume <- function(hh_size = 4, gpcd = 50){
vol = hh_size * gpcd * 30
return(vol)
}
calculate_bill <- function(volume = 6000){
fixed_charges = 7.63 + 9.85 + 1.8
vol_rate = 5.76 + 2.71
bill = fixed_charges+ (vol_rate*volume)/748.052
return(bill)
}
gradeAffordability <-function(HBI,PPI){
grade <- NA
grade <- case_when(
HBI < 7 & PPI < 20 ~ "Low Burden",
HBI < 7 & PPI>=20 & PPI <35 ~ "Moderate-Low Burden",
HBI < 7 & PPI >=35 ~ "Moderate-High Burden",
HBI >= 7 & HBI < 10 & PPI >= 35 ~ "High Burden",
HBI >= 7 & HBI < 10 & PPI >= 20 & PPI < 35 ~ "Moderate-High Burden",
HBI >= 7 & HBI < 10 & PPI < 20 ~ "Moderate-Low Burden",
HBI >= 10 & PPI >= 35 ~ "Very High Burden",
HBI >= 10 & PPI >= 20 & PPI < 35 ~ "High Burden",
HBI >= 10 & PPI < 20 ~ "Moderate-High Burden"
)
return(grade)
}
save.image("../cache/components_prep.RData")
load("../cache/components_prep.RData")
We adopt the approach in the figure below, as elaborated in this report.
HBI/PPI Matrix
#HBI Calculation
# Place
pl.w <- pl.w %>% mutate(bill_size_avg = calculate_bill(volume = calculate_volume(hh_size=hh_size_avg_overallE)))
pl.w <- pl.w %>% mutate(bill_size_1 = calculate_bill(volume = calculate_volume(hh_size=1)))
pl.w <- pl.w %>% mutate(bill_size_2 = calculate_bill(volume = calculate_volume(hh_size=2)))
pl.w <- pl.w %>% mutate(bill_size_3 = calculate_bill(volume = calculate_volume(hh_size=3)))
pl.w <- pl.w %>% mutate(bill_size_4 = calculate_bill(volume = calculate_volume(hh_size=4)))
pl.w <- pl.w %>% mutate(bill_size_5 = calculate_bill(volume = calculate_volume(hh_size=5)))
pl.w <- pl.w %>% mutate(bill_size_6 = calculate_bill(volume = calculate_volume(hh_size=6)))
pl.w <- pl.w %>% mutate(bill_size_7 = calculate_bill(volume = calculate_volume(hh_size=7)))
pl.w$HBI_size_avg <- 100*pl.w$bill_size_avg/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_1 <- 100*pl.w$bill_size_1/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_2 <- 100*pl.w$bill_size_2/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_3 <- 100*pl.w$bill_size_3/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_4 <- 100*pl.w$bill_size_4/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_5 <- 100*pl.w$bill_size_5/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_6 <- 100*pl.w$bill_size_6/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_7 <- 100*pl.w$bill_size_7/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$PPI <- 100*(pl.w$inc_poverty_countE-pl.w$inc_poverty_count_gr2E)/pl.w$inc_poverty_countE
pl.w$aff_grade_avg <- gradeAffordability(pl.w$HBI_size_avg,pl.w$PPI)
pl.w$aff_grade_1 <- gradeAffordability(pl.w$HBI_size_1,pl.w$PPI)
pl.w$aff_grade_2 <- gradeAffordability(pl.w$HBI_size_2,pl.w$PPI)
pl.w$aff_grade_3 <- gradeAffordability(pl.w$HBI_size_3,pl.w$PPI)
pl.w$aff_grade_4 <- gradeAffordability(pl.w$HBI_size_4,pl.w$PPI)
pl.w$aff_grade_5 <- gradeAffordability(pl.w$HBI_size_5,pl.w$PPI)
pl.w$aff_grade_6 <- gradeAffordability(pl.w$HBI_size_6,pl.w$PPI)
pl.w$aff_grade_7 <- gradeAffordability(pl.w$HBI_size_7,pl.w$PPI)
pl.bill<- pl.w %>% st_drop_geometry() %>% select(bill_size_avg,bill_size_1,bill_size_2,bill_size_3,bill_size_4,bill_size_5,bill_size_6,bill_size_7) %>%
pivot_longer(cols = c(bill_size_avg,bill_size_1,bill_size_2,bill_size_3,bill_size_4,bill_size_5,bill_size_6,bill_size_7),names_to="hh_size",names_prefix="bill_size_",values_to="bill")
pl.HBI<- pl.w %>% st_drop_geometry() %>% select(HBI_size_avg,HBI_size_1,HBI_size_2,HBI_size_3,HBI_size_4,HBI_size_5,HBI_size_6,HBI_size_7) %>%
pivot_longer(cols = c(HBI_size_avg,HBI_size_1,HBI_size_2,HBI_size_3,HBI_size_4,HBI_size_5,HBI_size_6,HBI_size_7),names_to="hh_size",names_prefix="HBI_size_",values_to="HBI")
pl.aff<- pl.w %>% st_drop_geometry() %>% select(PPI,hh_income_upper_limit_Q1E,aff_grade_avg,aff_grade_1,aff_grade_2,aff_grade_3,aff_grade_4,aff_grade_5,aff_grade_6,aff_grade_7) %>%
pivot_longer(cols = c(aff_grade_avg,aff_grade_1,aff_grade_2,aff_grade_3,aff_grade_4,aff_grade_5,aff_grade_6,aff_grade_7),names_to="hh_size",names_prefix="aff_grade_",values_to="aff_grade")
pl.aff <- pl.aff %>% left_join(pl.HBI,by="hh_size") %>% left_join(pl.bill,by="hh_size") %>% select(hh_size,bill,aff_grade,HBI,hh_income_upper_limit_Q1E,PPI)
print(paste0("The Average Household Size for Naperville is ",
round(pl.w$hh_size_avg_overallE,2)))
[1] “The Average Household Size for Naperville is 2.79”
print(paste0("At 50 GPCD, this corresponds to a monthly water volume of ",
calculate_volume(hh_size=pl.w$hh_size_avg_overallE)))
[1] “At 50 GPCD, this corresponds to a monthly water volume of 4185”
print(paste0("This corresponds to an water and sewer bill of $",
round(pl.w$bill_size_avg,2)))
[1] “This corresponds to an water and sewer bill of $66.67”
print(paste0("The Upper limit of the lowest quintile of annual household income for Naperville is $",
pl.w$hh_income_upper_limit_Q1E,
" or $",
round(pl.w$hh_income_upper_limit_Q1E/12,2),
" monthly"))
[1] “The Upper limit of the lowest quintile of annual household income for Naperville is $57958 or $4829.83 monthly”
print(paste0("The Poverty Prevalence Indicator (PPI) for Naperville overall is % ",
round(pl.w$PPI,1)))
[1] “The Poverty Prevalence Indicator (PPI) for Naperville overall is % 9.4”
print(paste0("By the HBI/PPI Matrix, the Water and Sewer Cost Burden in Naperville overall is graded below, assuming an average housheold size as well as household sizes 1-7"))
[1] “By the HBI/PPI Matrix, the Water and Sewer Cost Burden in Naperville overall is graded below, assuming an average housheold size as well as household sizes 1-7”
datatable(pl.aff) %>%
formatCurrency(c('bill','hh_income_upper_limit_Q1E'),'$') %>%
formatRound(c('HBI','PPI'),2)
Below we visualize HBI, PPI, and the overall affordability grade as calcualted at the tract level.
# Tract
# Place
tr.w <- tr.w %>% mutate(bill_size_avg = calculate_bill(volume = calculate_volume(hh_size=hh_size_avg_overallE)))
tr.w <- tr.w %>% mutate(bill_size_1 = calculate_bill(volume = calculate_volume(hh_size=1)))
tr.w <- tr.w %>% mutate(bill_size_2 = calculate_bill(volume = calculate_volume(hh_size=2)))
tr.w <- tr.w %>% mutate(bill_size_3 = calculate_bill(volume = calculate_volume(hh_size=3)))
tr.w <- tr.w %>% mutate(bill_size_4 = calculate_bill(volume = calculate_volume(hh_size=4)))
tr.w <- tr.w %>% mutate(bill_size_5 = calculate_bill(volume = calculate_volume(hh_size=5)))
tr.w <- tr.w %>% mutate(bill_size_6 = calculate_bill(volume = calculate_volume(hh_size=6)))
tr.w <- tr.w %>% mutate(bill_size_7 = calculate_bill(volume = calculate_volume(hh_size=7)))
tr.w$HBI_size_avg <- 100*tr.w$bill_size_avg/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_1 <- 100*tr.w$bill_size_1/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_2 <- 100*tr.w$bill_size_2/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_3 <- 100*tr.w$bill_size_3/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_4 <- 100*tr.w$bill_size_4/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_5 <- 100*tr.w$bill_size_5/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_6 <- 100*tr.w$bill_size_6/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_7 <- 100*tr.w$bill_size_7/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$PPI <- 100*(tr.w$inc_poverty_countE-tr.w$inc_poverty_count_gr2E)/tr.w$inc_poverty_countE
tr.w$aff_grade_avg <- gradeAffordability(tr.w$HBI_size_avg,tr.w$PPI)
tr.w$aff_grade_1 <- gradeAffordability(tr.w$HBI_size_1,tr.w$PPI)
tr.w$aff_grade_2 <- gradeAffordability(tr.w$HBI_size_2,tr.w$PPI)
tr.w$aff_grade_3 <- gradeAffordability(tr.w$HBI_size_3,tr.w$PPI)
tr.w$aff_grade_4 <- gradeAffordability(tr.w$HBI_size_4,tr.w$PPI)
tr.w$aff_grade_5 <- gradeAffordability(tr.w$HBI_size_5,tr.w$PPI)
tr.w$aff_grade_6 <- gradeAffordability(tr.w$HBI_size_6,tr.w$PPI)
tr.w$aff_grade_7 <- gradeAffordability(tr.w$HBI_size_7,tr.w$PPI)
hbi <- select(tr.w,
HBI_size_avg,
HBI_size_1,
HBI_size_2,
HBI_size_3,
HBI_size_4,
HBI_size_5,
HBI_size_6,
HBI_size_7,
)
ppi <- select(tr.w,
PPI)
aff <- select(tr.w,
aff_grade_avg,
aff_grade_1,
aff_grade_2,
aff_grade_3,
aff_grade_4,
aff_grade_5,
aff_grade_6,
aff_grade_7)
m1 <- mapview(aff,burst=TRUE) + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m2 <- mapview(hbi,burst=TRUE,col.regions=viridis::magma(7)) + mapview(b,alpha.regions=0,col.regions="green",stroke=TRUE,lwd=3,color="green",layer.name="Municipal Boundary")
m3 <- mapview(tr.w,zcol="hh_income_upper_limit_Q1E",col.regions=viridis::plasma(7),layer.name="Lower Quintile Income") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m4 <- mapview(ppi,burst=TRUE) + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
sync(m1,m2,m3,m4)
Below, we simulate cutoff rates as a probabilistic function of HBI and PPI, and Customer Assistance participation rates as a probabilistic function of cutoff rates and white population share. We simulate participation in conservation/ efficiency appliance incentive programs as a function of homeownership rates and rate of population with at least a bachelors degree.
In reality, one would hope to use actually utility records of participation in these programs at the address level, or aggregated to census tracts.
tr.w$white_perc <- tr.w$pop_race_whiteE/tr.w$pop_race_countE
P1 <- select(tr.w, NAME, GEOID, hh_type_total_countE, aff_grade_avg,
aff_grade_1,
aff_grade_2,
aff_grade_3,
aff_grade_4,
aff_grade_5,
aff_grade_6,
aff_grade_7,
HBI_size_avg,
HBI_size_1,
HBI_size_2,
HBI_size_3,
HBI_size_4,
HBI_size_5,
HBI_size_6,
HBI_size_7,PPI, hh_income_upper_limit_Q1E, hh_size_avg_overallE,white_perc, pop_owner_occupied_unitsE, pop_rental_unitsE,
pop_educ_attainment_countE, pop_educ_bachelorE, pop_educ_masterE, pop_educ_docE, pop_educ_profE)
set.seed(20212002)
#Cutoffs and CAP
z <- 5 -0.03* P1$HBI_size_avg - 0.02*P1$PPI - 0.04 * P1$HBI_size_avg*P1$PPI
P1$Cutoff_Perc = 100*(1/(1+exp(z)))
P1$CAP_Perc <- P1$Cutoff_Perc * (1-0.3*(1-P1$white_perc))
#Conservation (ever)
P1$bachAbove <- 100*(P1$pop_educ_bachelorE + P1$pop_educ_masterE + P1$pop_educ_docE + P1$pop_educ_profE)/P1$pop_educ_attainment_countE
P1$own_rate <- 100*P1$pop_owner_occupied_unitsE/(P1$pop_rental_unitsE+P1$pop_owner_occupied_unitsE)
z <- 4 - 0.01 * P1$bachAbove - 0.02 * P1$own_rate
P1$Incentive_rate <- 100*(1/(1+exp(z)))
m5 <- mapview(P1,zcol="Cutoff_Perc", layer.name="Account cutoff (%)") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m6 <- mapview(P1,zcol="CAP_Perc", layer.name="CAP participation (%)") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m7 <- mapview(P1,zcol="Incentive_rate", layer.name="Efficiency Incentive participation (%)") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
sync(m5,m6,m7)
We model all water quality events as basically random occurrences with probabilities between 0 and 10% and coverage of 50-100% at the tract level
set.seed(32896)
prob <- runif(length(P1$GEOID), min=0,max=0.1)
coverage <- runif(length(P1$GEOID),min=50, max=100)
event <- rbinom(length(P1$GEOID),1, prob=prob)
P1$Interrup <- event*coverage
set.seed(328435396)
prob <- runif(length(P1$GEOID), min=0,max=0.1)
coverage <- runif(length(P1$GEOID),min=50, max=100)
event <- rbinom(length(P1$GEOID),1, prob=prob)
P1$boil_rate <- event*coverage
This would be a Utility-Level measure only, and is available from EPA SDWIS
In this case, the API call for Naperville involves the PWSID IL0434670:
https://enviro.edap-cluster.com/cache/query/enviro3/efservice/VIOLATION/PWSID/IL0434670/IS_HEALTH_BASED_IND/Y/ROWS/0:100/JSON
In this case, there are no health-based violations in the past 10 years, so the value of the indicator is 0.
Below we export the Pillar 1 indicators.
P1_tract <- select(P1,-pop_owner_occupied_unitsE,-pop_educ_attainment_countE,-pop_rental_unitsE,-pop_educ_bachelorE,-pop_educ_masterE,-pop_educ_docE,-pop_educ_profE)
st_write(P1_tract,"../out/P1_tract.geojson",append=FALSE)
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Error 6: DeleteLayer() not supported by
## this dataset.
## Deleting layer not supported by driver `GeoJSON'
## Deleting layer `P1_tract' failed
## Writing layer `P1_tract' to data source `../out/P1_tract.geojson' using driver `GeoJSON'
## Updating existing layer P1_tract
## Writing 45 features with 30 fields and geometry type Multi Polygon.
write.csv(st_drop_geometry(P1_tract),"../out/P1_tract.csv")
P1_overall <- P1_tract %>% st_drop_geometry() %>% summarise(across(c(HBI_size_avg:hh_income_upper_limit_Q1E,Incentive_rate:boil_rate,Cutoff_Perc:CAP_Perc),~weighted.mean(.,w=hh_type_total_countE))) %>% mutate(health_based_sdwa_violations=0)
write.csv(P1_overall,"../out/P1_overall.csv")
P2_tract <- tr.w %>%
select(NAME,GEOID,hh_type_total_countE) %>%
mutate(staff_count = round(runif(length(NAME),min=0,max=10),0),
staff_per_1000hh = 1000*staff_count/hh_type_total_countE)
P2_overall <- P2_tract %>% st_drop_geometry %>% summarize(staff_count_in_service_area = sum(staff_count),
staff_per_1000hh = sum(staff_count)/sum(hh_type_total_countE))
P2_overall$staff_total<-300
Diversity requirements in hiring
Diversity requirements in procurement
Local requirements for procurement/contracting
Participation in local education, training and workforce development programs
set.seed(35455)
P2_overall <- P2_overall %>%
mutate(require_diverse_hiring = rbinom(1,1,0.5),
require_diverse_procurement = rbinom(1,1,0.5),
require_local_procurement = rbinom(1,1,0.5),
regional_training_education_programs = round(runif(1, min=0,max=5)),
region_training_education_participation = round(runif(1,min=0,max=100))
)
Below we export the Pillar 2 inidicators in a standardized format.
st_write(P2_tract,"../out/P2_tract.geojson",append=FALSE)
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Error 6: DeleteLayer() not supported by
## this dataset.
## Deleting layer not supported by driver `GeoJSON'
## Deleting layer `P2_tract' failed
## Writing layer `P2_tract' to data source `../out/P2_tract.geojson' using driver `GeoJSON'
## Updating existing layer P2_tract
## Writing 45 features with 5 fields and geometry type Multi Polygon.
P2_tract <- st_drop_geometry(P2_tract)
write.csv(P2_tract,"../out/P2_tract.csv")
write.csv(P2_overall,"../out/P2_overall.csv")
Number Community water planning meetings (in each Tract)
Proportion structures with lead service lines (by tract)
Water main breaks/length pipe (by tract)
Sewer overflows (by tract)
Water quality complaints received/ addressed (by tract)
Lead service line replacements (proportion of lead service lines by tract)
Number leaks detected (by tract)
Expenditure incurred for detection measures
Infrastructure repair/rehabilitation/replacement expenditure (by tract)
Real losses reduced due to repairs/rehab/replacement (by tract)
P3 <- tr.w %>%
select(NAME,GEOID,hh_type_total_countE)
P3.c <- st_join(P3,lead.c)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
P3.c <- P3.c %>% mutate(water_quality_complaints_received = rbinom(length(GEOID),1,0.08),
water_quality_complaints_resolved = water_quality_complaints_received*rbinom(length(GEOID),1,0.9),
ft_mains = runif(length(GEOID),55,164),
lead_line = as.numeric(SUBTYPE=="LEAD"))
P3_tract <- P3.c %>% st_drop_geometry() %>%
group_by(GEOID,NAME,hh_type_total_countE) %>% add_tally() %>%
summarize(conn = mean(n),
water_quality_complaints_received = sum(water_quality_complaints_received),
water_quality_complaints_resolved = sum(water_quality_complaints_resolved),
ft_mains = sum(ft_mains),
lead_lines = sum(lead_line,na.rm=TRUE)) %>%ungroup()
## `summarise()` has grouped output by 'GEOID', 'NAME'. You can override using the `.groups` argument.
P3_tract <- P3_tract %>%
mutate(meeting_count = round(runif(1, min=0,max=4)),
meeting_count_1000hh = meeting_count/hh_type_total_countE,
main_breaks = floor(runif(length(GEOID),min=0,max=3)),
main_breaks_1000ft_pipe = main_breaks/(ft_mains/1000),
prop_lead_lines = lead_lines/conn,
leaks_detected = floor(runif(length(GEOID),min=main_breaks, max = 10)),
prop_lead_lines_replaced = case_when(lead_lines > 0 ~ runif(1,min=0,max=1),
lead_lines == 0 ~ NA_real_),
sewer_overflows = floor(runif(length(GEOID), min=0, max=2)),
leak_detection_expenditure = runif(length(GEOID),min=0,max=1000000),
network_renewal_percent = runif(length(GEOID),min=0,max=20),
network_renewal_expenditure =runif(length(GEOID),min=0,max=1000) * network_renewal_percent,
real_losses_percent_delivered = runif(length(GEOID),min=0,max=100),
real_loss_reduction_percent = case_when(real_losses_percent_delivered > 5 ~ runif(length(GEOID),min=0,max=95),
real_losses_percent_delivered <= 5 ~ runif(length(GEOID),min=0,max=50))
)
Utility Resilience Index (URI)
Infrastructure Leakage Index (by DMA -> summarize over Census Tract if possible)
P3_overall <- P3_tract %>%
summarize(hh = sum(hh_type_total_countE),
conn = sum(conn),
meeting_count = sum(meeting_count),
main_breaks = sum(main_breaks),
ft_mains = sum(ft_mains),
lead_lines = sum(lead_lines),
leaks_detected = sum(leaks_detected),
lead_lines_replaced = mean(prop_lead_lines_replaced,na.rm=TRUE),
sewer_overflows = sum(sewer_overflows)
# network_renewal_percent2 = weighted.mean(x=network_renewal_percent,w=ft_mains,na.rm=TRUE),
#real_losses_percent_delivered2 = weighted.mean(x=real_losses_percent_delivered,w=conn,na.rm=TRUE),
#real_loss_reduction_percent2 = weighted.mean(x=real_losses_reduction_percent,w=conn,na.rm=TRUE)
)
P3_overall <- P3_overall %>% mutate(
URI = floor(runif(1, min=0, max=100)),
ILI = runif(1,min=1,max=4),
prop_lead_lines = lead_lines/conn,
breaks_per_1000ft = main_breaks/(ft_mains/1000),
leaks_per_1000ft = leaks_detected/(ft_mains/1000),
sewer_overflows_per_1000conn = sewer_overflows/conn
)
P3_overall$network_renewal_percent = weighted.mean(x=P3_tract$network_renewal_percent,w=P3_tract$ft_mains,na.rm=TRUE)
P3_overall$real_losses_percent_delivered = weighted.mean(x=P3_tract$real_losses_percent_delivered,w=P3_tract$conn,na.rm=TRUE)
P3_overall$real_loss_reduction_percent = weighted.mean(x=P3_tract$real_loss_reduction_percent,w=P3_tract$conn,na.rm=TRUE)
P3_tract <- left_join(P3,P3_tract,by=c("GEOID","NAME","hh_type_total_countE"))
Below we export the Pillar 3 indicators
st_write(P3_tract,"../out/P3_tract.geojson",append=FALSE)
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Error 6: DeleteLayer() not supported by
## this dataset.
## Deleting layer not supported by driver `GeoJSON'
## Deleting layer `P3_tract' failed
## Writing layer `P3_tract' to data source `../out/P3_tract.geojson' using driver `GeoJSON'
## Updating existing layer P3_tract
## Writing 45 features with 21 fields and geometry type Multi Polygon.
P3_tract <- P3_tract %>% st_drop_geometry()
write.csv(P3_tract,"../out/P3_tract.csv")
write.csv(P3_overall,"../out/P3_overall.csv")
Below you can compare any two attributes side-by-side by selecting the appropriate layers in the layer control icon of each map in the upper left.
tracts <- left_join(P1_tract,P2_tract,by=c("GEOID","NAME"))
tracts <- left_join(tracts,P3_tract,by=c("GEOID","NAME"))
tracts <- select(tracts,-GEOID,-NAME)
t1 <- mapview(tracts,burst=TRUE,homebutton=FALSE) + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
sync(t1,t1,ncol=1)
save.image("../cache/indicators_ungraded.RData")